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Abstract. We formulate an energy-based atomistic- to-continuum coupling 
method based on blending the quasicontinuum method for the simulation of 
crystal defects. We utilize theoretical results from [24 32 to derive optimal 
choices of approximation parameters (blending function and finite element 
grid) for microcrack and di-vacancy test problems and confirm our analytical 
predictions in numerical tests. 

1. Introduction 

A major goal of materials science is to predict the macroscopic properties of ma- 
terials from their microscopic structure. For this purpose, it is necessary to under- 
stand the behavior of defects in these materials. We propose a computational tool, 
the energy-based blended quasicontinuum method (BQCE), for simulating defects 
such as cracks, dislocations, vacancies, and interstitials in crystalline materials. 

Accurate modeling of the region near a defect requires the use of computationally 
expensive atomistic models. Such models are practical only for small problems. 
However, a defect may interact with a large region of the material through long- 
range elastic fields. Thus, accurate simulation of defects requires the use of a large 
computational domain; typically, the size required rules out the use of atomistic 
models for the entire region of interest. 

Fortunately, the long-range elastic fields generated by a defect are well described 
by continuum models which can be efficiently computed using the finite-element 
method. Thus, defects can be accurately and efficiently simulated by coupled mod- 
els which use an atomistic model near the defect and a continuum model elsewhere. 
We call any such model an atomistic-to- continuum coupling. 

Many atomistic-to-continuum couplings have been proposed in recent years [I]- 
[6|[l3j[T6j[2lJ|28|[30]; see 19 31 for a survey of atomistic-to-continuum couplings 



and computational benchmark tests. These couplings fall into two major classes: 
energy-based and force-based. Energy-based couplings provide an approximation 
to the atomistic energy of a configuration of atoms. Force-based couplings provide a 
non-conservative force-field which approximates the forces on each atom under the 
atomistic model. Our BQCE method is an energy-based coupling. Both types of 
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couplings have intrinsic advantages; the development of energy-based couplings is 
especially important for finitc-tempcrature applications since equilibrium statistical 
properties and transition rates can be directly approximated 10 18 . 

Here, and throughout, we only consider concurrent coupling methods. An al- 
ternative approach are upscaling methods such as 26 , which achieve increased 
resolution near defects by adding finer scales to a coarse scale description rather 
than by decomposing space into fine (atomistic) and coarse (continuum) scale mod- 
els. 

The primary source of error for most (concurrent) energy-based couplings is the 
ghost force. We say that a coupling suffers from ghost forces if it predicts non- 
zero forces on the atoms in a perfect lattice. Although many attempts have been 
made to develop an energy-based coupling free from ghost forces such couplings are 



currently known only for a limited range of problems 11 25 28 30 
Shapeev's method 



28 



applies to one and two-dimensional simple crystals with 
an atomistic energy based on a pair interaction model, and can be extended to 3D 
if a modified "continuum model" is used 29 . GR-ACQwas proposed in TlpO , and 
has recently been implemented for a two dimensional crystal with nearest neighbor 
many-body interactions in [25] . No ghost force free methods are currently known 
for three-dimensional crystals, for multi-lattice crystals (except in ID 22 ), or for 



atomistic models with general many-body interactions. The field-based coupling of 
Iyer and Gavini |16| is another interesting approach; however, it is unclear at present 
whether it is competitive in terms of computational complexity and generality. 

In our BQCE method, the ghost forces cannot be eliminated but can be con- 
trolled in terms of an additional approximation parameter (the blending width) [24| 
32 . BQCE applies to a wide range of problems for which no ghost force free methods 
are known; these problems include three-dimensional crystals with general many- 
body interactions as well as multi-lattices. This makes it an attractive method for 
such challenging and physically important problems. 

The key feature of BQCE is a blending region where the atomistic and continuum 
contributions to the total energy are smoothly mixed. The ghost forces of the 
BQCE method can be made arbitrarily small by increasing the size of this blending 
BQCE shares the idea of a blending region with the bridging domain 



24,32 



region 

method [6], the AtC coupling [I], and the Arlequin method |5j. By contrast, the 
energy-based quasicontinuum method (QCE) [21] , Shapeev's method [28] , and GR- 
AC 11 25 30 exhibit an abrupt transition between the atomistic and continuum 



models. We call any method with a blending region a blended method, and we call 
the weights which mix the atomistic and continuum contributions to the energy a 
blending function. 

Both the bridging domain method and AtC coupling are very general formula- 
tions, each of them incorporating BQCE and QCE as special cases. Our BQCE for- 
mulation provides a set of specific instructions for the successful implementation of 
a blended method. We identify two important practical differences between BQCE 
and the bridging domain and AtC coupling methods. First, the BQCE method 
specifies a strong coupling between the atomistic and continuum regions, whereas 
weak couplings based on Lagrange multipliers or the penalty method have been used 
in most work involving the bridging domain and AtC coupling methods. Second, in 



lr The acronym "GR-AC" was introduced in 
atomistic-to-continuum coupling method." 



251. It stands for "geometry reconstruction-based 
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BQCE, we blend the atomistic site-energy with a continuum site-energy based on 
the continuum site-energy defined in some formulations of the QCE method (see 
Section [3]). This guarantees that BQCE correctly predicts the total energy of a 
perfect lattice subjected to uniform strain. 



Our approach to blending is supported by rigorous analysis. In 32 , we showed 
that the ghost force error of BQCE in ID does indeed decrease with the size of 
the blending region, and we found that the error is minimized when the blending 
function is a Hermite cubic-spline. We also found that the error of the BQCE 
method in predicting lattice instabilities can be reduced by increasing the size of 
the blending region. 



In this paper, we have extended the atom-based formulation (4.2 1 given for 



ID problems in [32] to a volume-based multiD formulation (4.3) that allows finite 
element coarse-graining in the continuum region. 

The implementation of BQCE requires the choice of two approximation param- 
eters: a blending function [3 and a finite-element mesh T which is used to compute 
the continuum contribution to the energy. In Section |4.4[ we give optimal choices 
of (3 and T to minimize global error norms for the problem of a point defect in a 
2D crystal, based on theoretical results in [24] . 

We demonstrate the validity of these results in computational test problems in 
which we simulate a microcrack (see Figure [T]) and a di- vacancy. 

Outline. In Sections [2] and [3] we introduce an atomistic model for simple lattices 
with general many-body interactions and its corresponding continuum Cauchy- 
Born approximation. In Section [4] we give a precise formulation of BQCE in this 



setting. In Section 4.4 we offer guiding principles on choosing the blending function 



(3 and the mesh T based on our analysis in 24 . The results of this analysis are 
summarized in Table [l] In Section [5] we describe the details of our numerical 
experiment. We use an atomistic energy based on the Embedded Atom Method. 
We did not choose our atomistic energy to model any specific physical material; 
instead, the atomistic energy is a toy model chosen for its simplicity. The observed 



rates of convergence are in agreement with the rates predicted in Section 4.4 In 

1 9 — - — 

particular, the error of BQCE in the W • semi-norm decreases as DoF 2 where 

DoF is the number of degrees of freedom. We summarize our results and discuss 
open questions concerning local quantities of interest in the concluding Section [6] 

2. The Atomistic Energy 

Let A be a e?-dimensional Bravais lattice. We call A the reference lattice. In the 
present work, we consider only monatomic crystals. That is, we assume that each 
site of the lattice A is the reference position of a single atom, and that all atoms 
arc identical. Let 

9 := {y : A -> K rf ; y(0 f y( V ) for all £ f 77} 

be the set of deformations of A. Let f2 a C A be a finite subset of A. We call fl a 
the atomistic computational domain. 

In the Embedded Atom Method (EAM) [§], the energy of O a subjected to de- 
formation y takes the form 

(2.1) E\y) := E { E ~ ViOl) + ^( £ p(\vW ~ V®\)) }> 
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where tj> is a pair potential, p is an electron density function, and G is an embedding 
function. We call the inner sum 

(2.2) q(y) := £ U(\y{ V ) - y(OI) + g( £ - 2/(01) 



2 

the atomistic site-energy of atom £. We observe that the sum defining the energy of 
an atom is finite in practice even though summation ranges over the infinite lattice 
A. This is because the pair potential <p an d electron density function p are taken 
to have a cut-off radius r c such that 

4>{r) = p(r) — for all r > r c . 

For defining the atomistic model and the BQCE method, we can in fact consider 
a more general class of potentials than EAM. We only require that the total energy 
can be decomposed into a sum of localized site-energies associated with each atom. 
By localized, we mean that the energy associated with an atom £ does not depend 
on the positions of atoms beyond a certain cut-off distance. Such an assumption 
may be violated for certain energies arising from quantum mechanics, but does 
hold for most empirical potentials including EAM, the bond-angle potentials, and 
so forth. In addition, we require that the site-energies are homogeneous; that is, 
the energies of atoms which have the same local environment are the same. 

To make these assumptions precise, we let £f(y) denote the site-energy associated 
with atom £ under deformation y and require that it is of the form 

(2.3) £f(y) := V({y( V ) - y(® : ry e A \ {£}, \y(v) ~ y(0\ < r c }), 

where V is the site-energy potential. We assume that the resulting site-energies are 
twice continuously differentiable, that is, £f £ C 2 {f3f\ The restricted dependence 
of £? on atoms within the cut-off radius may also be expressed in the form 

d£f 

(y) = for all i) with 1 3/(77) - y(£)| > r c . 



dy{ri) 

This quantifies the requirement that the site-energy is localized. 

Given £^ we define the energy of f2 a subjected to a deformation y by 

(2-4) £ a (2/) := ]T q(y). 



We call an energy of the form (2.4) where f| satisfies (2.3 1 homogeneous. When we 
define the Cauchy-Born strain energy density corresponding to (2.4) in Section |3j 
we use that £ a is homogeneous: if £ a is not homogeneous, the energy per unit 
volume in a perfect lattice subjected to uniform strain may be more difficult to 
obtain (see [l] and references therein for examples). 

Remark 2.1. The locality assumption can, in principle, be replaced by an assump- 
tion that the interaction strength decays sufficiently rapidly with increasing distance 
between atoms. 

Remark 2.2. Homogeneity of the site-energy is our main assumption that is vio- 
lated for multi-lattices. We show in \2J^ how to generalize our formulation for that 
scenario. 
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3. The Cauchy-Born Site Energy 



BQCE is a coupling of an atomistic energy based on (2.4 1 with a coarse-grained 
continuum elastic energy based on the Cauchy-Born strain energy density (3.1). 
Let vor(£) denote the Voronoi cell of a site £ G A. Then the Cauchy-Born strain 
energy density W : M. dxd — > K U {±00} corresponding to V is defined by 

(3-D W(F) := jJ^sA 

where y F G is the homogeneous deformation y (£) — F£. W(F) may be inter- 
preted as the energy per unit volume in A subjected to strain F. Note that the 



assumption of homogeneity (2.3) ensures that £f{y ) — £§(y ) for all £ G A 



We will use the Cauchy-Born strain energy density (3.1) to derive a coarse- 
grained continuum energy suitable for coupling with the atomistic energy (2.4). 
First, we define a space of coarse-grained deformations. Let A rep C A denote a 
set of representative atoms (repatoms), let T be a triangulation with vertices A rep , 
and let P 1 (T) denote the set of all functions y : M. d — » M. d that are continous 
and piecewise affine with respect to T. We call P 1 (T) the set of coarse-grained 
deformations. The Cauchy-Born energy of a deformation y G P 1 (T) in a domain 
Q is then given by 

£ c (y) := / W{Vy{x))dx. 



<> 



The Cauchy-Born approximation is analyzed, for example, in |7l[T2 15 



The definition of the QCE method 21 and our construction of the BQCE 



method in the next section use a Cauchy-Born site-energy £?, which is analogous 
to the atomistic site-energy For y G -P 1 (T) and £ G A, we define £| by 



£l(y)--= / W{Vy{x))dx 
Jvor(£) 

(3.2) = |vor(Onr|W(Vy|T). 



In formula (3.2), |vor(£) D T\ denotes the volume of the intersection of vor(£) with 



the element T. We observe that the sum on the right hand side of equation (3.2) 
is finite because only finitely many elements T G T can intersect vor(£). 

4. The Blended Quasicontinuum Energy 

4.1. Formulation of the BQCE method. The BQCE method is an atomistic- 
to-continuum coupling based on the QCE method of Tadmor et al. 21 . In QCE, 
the reference domain f2 a is partitioned into an atomistic region A and a continuum 
region C, and the QC energy £ qc : P 1 (T) — > K is defined by 

(4.1) ^j-Ew+Ew 

In BQCE, the atomistic and continuum energies per atom are weighted averages. 
Given a blending function /3 : O a -> [0,1] the BQCE energy £ p : P^T) -> K is 
defined by 

(4-2) £ (y) := X! ^(6^ c (y) + (1 - /3(0)^ a (2/). 
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We observe that the QCE energy with continuum region C is the same as the BQCE 
energy with (3 chosen as the characteristic function of C. Our formulation of BQCE 
is similar in spirit to the bridging domain method [6], the AtC coupling El, and 
the Arlequin method (5l; it differs from [4H6] in that we specify a strong coupling 
between the atomistic and continuum models, and accordingly our method is based 
on the minimization of an energy. Moreover, our formulation in terms of atomistic 
and continuum energies per atom guarantees that the BQCE energy is exact under 
homogenous deformations. 

The BQCE energy can be rewritten in the form 

e p (v) = E #0 i w{ yy) dx+(i- mmiy) 

-/vor(f) 

= E E /5(0|vor(0 R T\W(Vy\ T ) + E O 1 ~ 

(4.3) = e 4w(yy\ T ) + E (! - ftt)) £ t(v)> 

TeT ten* 
where the BQCE- effective volume of the element T is defined by 

(4.4) 4:= E /3(£)|vor(£)nT|. 

Remark 4.1. The triangulation T need not cover the entire domain f2 a . For T a 
triangulation which covers only part ofR d , define 

Q(T) : = UrerT, and 

pl 0~) : = {y ■■ O(T) U A ->• R d : y piecewise affine w.r.t. T on fl(T)}. 

We observe that £^{y) is defined for y £ P X (T) if for every £ € fi a smc/i i/iai 
> we /lave vor(^) C 0(7^. /n particular, it is not necessary to assume that 
the triangulation T is refined to atomistic scale anywhere in the domain f2 a . It is 
possible that the use of a mesh which is not refined to atomistic scale may make the 
implementation of BQCE easier and more efficient in some cases. 



Remark 4.2 (Multi-lattices). If one interprets a multi-lattice as a simple lattice 
with a basis (see, e.g., f§ 24 ~3^), then the formulation of the atomistic energy (2.4) 
and of the BQCE energy (4.2) require no changes. The only difference is that in 
the general multi-lattice case, the site energy £"'(•'/) has fewer symmetries (however, 
see I '3 '3y for examples with high symmetry). 



4.2. Far-field boundary conditions. A typical application of the BQCE method 
is the simulation of a defect or defect region in an infinite crystal. To that end, 
we require far-field boundary conditions at the domain boundary. We propose two 
choices. 

4.2.1. Dirichlet boundary conditions. Let S C T be a finite subset of T, and let 
VL(S) := UsesS be a polygonal domain. When Dirichlet boundary conditions are 
imposed, the deformation of the boundary dfl(S) of f2(<S) is fixed to agree with 
some yo E P 1 (T) ■ Precisely, we let 

(4.5) Adm := {y G P^T) : y(x) = y (x) for all x £ dQ(S)} 
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denote the space of admissible deformations. We then solve the problem 

(4.6) Find y G argminf 13 (z), 

zGAdm 

where we interpret argmin zgAdm S^{z) as the set of local minimizers of £^ . 

4.2.2. Periodic boundary conditions. A popular method to construct artificial far- 
field boundary conditions is to formulate the problem in a periodic cell. To that end, 
suppose that f2(T) = T xs periodic, and let S C T be the finite element mesh 
on one periodic cell. That is, suppose that for some nonsingular matrix A G M. dxd 
whose columns are elements of A, we have 

T= (J {An + S}, 

and that the union is disjoint. 

Let F € WL dxd be nonsingular, and call F a macroscopic strain. Given F, we 
define the admissible set as 

(4.7) Adm := {y G P 1 (T) : y(x + An) = y(x) + F An for all x E 

The associated variational problem can again be stated as ( |4.6[ ). 

The columns of the matrix A should be understood as the directions in which the 
displacement u{x) = y(x) — Fx is periodic. Equivalently, the set {Ax : x G [Q, 
is a periodic cell. The macroscopic strain F should be understood as the average 
strain imposed on the crystal. If the size of the periodic cell is large, periodic 
boundary conditions with macroscopic strain F approximate far-field boundary 
conditions with uniform strain F imposed at infinity. 

Note that, while we only discuss a static formulation, we see no obstacle that 
prevents both our formulation and analysis to apply also to quasi-static problems, in 
which the energy is minimized at each loading step. However, we stress that further 
careful analysis is required to understand the accuracy, in particular stability, of 
the BQCE method as the deformation approaches bifurcation points (e.g., as in the 
formation or motion of crystal defects or the propagation of a crack) . 



4.3. Degrees of freedom (DoF). In Section 4.4 we present error estimates for 
the BQCE method in terms of the computational complexity, by which we mean 
the computational cost to compute the BQCE energy and its gradient. This is 
proportional to the number of degrees of freedom of the space of admissible functions 
Adm. Similarly, in Section [5] we plot errors against the number of degrees of 
freedom in the computed solution. 

The number of degrees of freedom (DoF) is simply the size of the nonlinear 
system that we solve to minimize the BQCE energy. In terms of the notation 
introduced in the previous section, this is the number of finite clement nodes in the 
triangulation S plus the number of lattice sites in A \ f2(<S) for which /3(£) < 1. 

In the reduced atomistic computation, which we describe in Section [4.5. 2| DoF 
denotes the number of unconstrained lattice sites. 



4.4. Computational complexity and optimal parameters. In 24 , we con- 
jecture an error estimate for BQCE in 2D and 3D. The conjecture is based on our 
error analysis of BQCE in ID [32) and on the consistency estimates for BQCE in 



2D and 3D in 24 . Following 23 Sec. 7.1], we use our conjectured error estimate 



to derive complexity estimates, which are bounds on the error of the BQCE method 
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in terms of the number of degrees of freedom. We use our complexity estimates to 
guide the choice of optimal approximation parameters T and /3. 

Let y a be a local minimizer of the atomistic energy, and let yp be a local minimizer 
of BQCE with the same boundary conditions as y a . Let h be the mesh size function 
of the triangulation T (h(x) — diam(T) for a.e. x € T), and let (3 be the blending 
function. Then we conjecture an estimate of the form: 

, , Err := ||Vy a - Vyp\\z, < \\hV 2 y a \\ LP(c) + ||V 2 /?|| L p 

( ' ' =: CG + GF, 

for given p e [l,oo]. For p — 2, we have proven this result in [24] ; for p ^ 2, it 
would be difficult to establish this result rigorously; and for p £ {1, oo}, it is unclear 



in what generality the result holds. In (4.8), V y a and V /3 should be interpreted 



as the second derivatives of smooth interpolants of y a and /3, and C :— supp(/3). 
The first term, CG = \\h'V 2 y ll \\Lp(c)> is the finite element coarsening error, while 
the second term, GF = ||V 2 /3||lp, measures the effect of the ghost forces. 

We now describe the problem of a point defect in a 2D crystal. To quantify the 
notion of a point defect, we assume that for some a > 0, 

(4.9) |VVO)l - r~ a , where r = \x\. 

It has been observed in numerical experiments that a — 2 for a dislocation, and 



that a = 3 for a vacancy 14 23 



We assume that the reference domain f2 is a roughly circular region of radius 
N atomic spacings centered at the origin. We let Kq > be the radius of the 
atomistic region surrounding the defect, and we let Ki > be the width of the 
blending region. We then choose a radial blending function /? of the form 

f0 ifr<X , 

(4.10) P(x) := I 0o (^ip) XK <r<K + Ki, 

[ 1 if K + Ki < r, 

where (3q : [0, 1] — > [0, 1] is a twice continuously differentiable function with /3q(0) = 
#(0) = P' (l) = and = 1. 

We summarize our complexity estimates for BQCE in Table [T] These estimates 



are proved in |24|. We distinguish three cases based on the value of 7 := In the 



second column, we give the optimal rates of convergence for BQCE. The optimal 
rates are attained when K\ is given in terms of Kq by the formula appearing in the 
third column and when the mesh size function h is given by 



(4.11) h(x) 



K 



Remark 4.3. When 7 > 1 and a > 2, another choice of parameters results in 
convergence at the optimal rate but with a better constant: we use 

a - 1 

K l= K% for A*=2"-f 

v 

in our numerical experiments below. This choice guarantees that the finite element 
error decreases at the same rate as the ghost force error. See '24-1 for additional 
suggestions related to choosing parameters. 
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Table 1. Complexity estimates and optimal approximation pa- 
rameters for BQCE 24 . The estimates above are for the problem 
of a point defect in a two-dimensional crystal. The variables K , 

ap 



4.4 



and 7 



p+2 ■ 



In all cases, 

7 



Ki, N, and a are defined in Section 

the optimal rate of convergence is attained when h = (jg^) an d 
when Ki is given in terms of K by the formula in the third column 
above. 



Case 


Complexity Estimate 




Optimal Parameters 


7 > 1 


||Vy a -V^|| iP <DoF m ^i- 


-i^-f} 


Kt=K 


7 = 1 


UVya-V^IU* <DoF max U" 


-h-h} 




7 < 1 


l|Vy a -V y/3 || L p<DoF raax a- 


-w> 


K l = K^N 1-1 



Remark 4.4. To make the explicit calculation of optimal parameters possible, we 



have made simplifying assumptions on the forms of (3 and h. In (4.10), we assume 



that (3 is radial. We explain how to construct (3 given general atomistic, continuum, 



and blending regions in Section 5.2 below. Similarly, we assume that h is continuous 
even though a nonconstant continuous function cannot be the mesh size function 
of a triangulation. In practice, we implement BQCE using a triangulation whose 



mesh size function approximates the continuous h defined in (4.11); we explain how 
to do this in Section \5.2\ We expect that the complexity estimates in Table [7] will 
hold for practical choices of (3 and h, not only for the simplified f3 and h used in 
the derivation of the estimates. 

Remark 4.5. The rates of convergence depend on the geometry of the problem 
and on the norm in which the error is measured. We observe that the error of 
BQCE does not decrease with DoF when measured in the W ' -seminorm; see the 
case p = 1 in Table^ This is because the W -1 ' 1 -norm of the ghost force does not 
decrease as the size of the blending region increases. We refer the reader to \2J$ for 
further discussion of the convergence of BQCE for other geometries and in other 
norms. 

Remark 4.6. We have proved that linear blending functions have a reduced rate 



of convergence compared to smooth blending functions of the form (4.10) (see \32\ 
for ID and ^ for 2D). In Table [§ we give comparative convergence rates for 
the problem of a point defect in a two-dimensional crystal with decay a = 3. We 
emphasize, however, that this is only an asymptotic estimate. Indeed in our micro- 
crack experiment in Section \5. 3\ we observe a significant preasymptotic regime where 
the two blending functions yield comparable errors. 

Remark 4.7. Although our subsequent numerical experiments are restricted to de- 
fects with zero Burgers vector and consequently 7 > 1, we consider more general 
situations in the above analysis, in order to emphasize the generality of our ap- 
proach. The implementation of a/c couplings for defects with non-zero Burgers 
vector poses additional challenges (what should one use for the reference domain in 
the neighborhood of a single dislocation?) that we will address in future work. As 
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linear BQCE 

||Vl/a- V^||lp 

< DoF"2 

< DoF"5 



smooth BQCE 

liV^a - II LP 

< DoF^s 

< DoF -1 



Table 2. Complexity estimates comparing the optimized linear 
and smooth BQCE method for the problem of a point defect in a 
two-dimensional crystal with decay a = 3. 



a matter of fact, we are not aware of numerical results from a/c coupling methods, 
simulating a single dislocation core in an otherwise defect-free lattice (as opposed 
to the simulation of a dislocation pair with total zero Burgers vector). 

4.5. Comparison of BQCE with other methods. We use the estimates above 
to compare the complexity of BQCE with direct atomistic simulation and with 
patch test consistent couplings, such as Shapeev's method. 



4.5.1. Shapeev's method and other patch test consistent methods. When a is small 
(slow decay of the elastic field), BQCE converges at the same rate as patch test 
consistent methods such as Shapeev's method. In fact, if a < 2 (e.g., a dislocation), 
then the VF^-error of Shapeev's method decreases as DoF~ 3 with the number of 
degrees of freedom [23]. We predict the same rate of convergence for BQCE when 
a < 2. On the other hand, when a is larger, patch test consistent couplings may 
converge faster than BQCE. When a > 2 (e.g., a vacancy, microcrack, or dislocation 
dipole), the W 1,2 -error of Shapeev's method decreases as DoF 5_T , whereas the 
W 1,2 -error of BQCE decreases as DoF~ 3 . Roughly speaking, this is due to the 
fact that, for small a, the coarse-graining error dominates, whereas for large a, the 
ghost force error dominates. 

According to the estimates above, Shapeev's method always converges at least as 
fast as BQCE. However, we remind the reader that patch test consistent methods 



are currently known only for a limited range of problems 11 25 28 30 . As far 



as we are aware, BQCE is the only convergent energy-based method for 2D or 
3D crystals with general many-body interactions. Moreover, even when patch test 
consistent methods are known, BQCE may be an attractive method for slowly 
decaying defects, since it converges at the same rate as patch test consistent methods 
for these problems. See [24] for a more detailed analytic comparison of BQCE with 
Shapeev's method. 



4.5.2. Reduced atomistic computations (ATM). In our numerical experiments be- 
low, we compare the performance of BQCE with the reduced atomistic computation 
(ATM), in which far-field boundary conditions are simulated by imposing Dirichlet 
boundary conditions directly on the atomistic model. Precisely, let 57 C A be an 
atomistic computational domain, and let F be a macroscopic strain as in Section [4.2| 
We then define the admissible space 



Adm := {y : A -> R d : = F£ for all £ € A \ SI}, 
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BQCE 


ATM 


2 < a 


ATM rate vs. 




\\Vy a - VvpWl* 


l|Vy a -Vy^|| L p 




BQCE rate 








a > 3 


higher 


p = 2 


< DoF"5 


< DoF _ 5 + 5( 3 ~ Q ) 


a = 3 


same 








a < 3 


lower 








a > 3 


higher 


p = oo 


< DoF" 1 


<DoF- 1+ ^ 3 - Q ) 


a = 3 


same 








a < 3 


lower 



Table 3. Complexity estimates for the optimized BQCE and 
ATM for the problem of a point defect in a two-dimensional crystal 
with decay a > 2 



24 . The variable a is defined in Section 4.4 



and we solve the variational problem 

(4.12) Find y £ argmin£ a (z). 

zGAdm 

We regard ATM as the simplest means of simulating far-field boundary condi- 
tions, and we feel that an atomistic-to-continuum coupling is successful only when 



it is more accurate than ATM. In 24 , we show that 



\Vy% - Vy,\\ LP < Nt +1 - a ~ DoFF + 5-f 



where N is the diameter of f2, y^ solves (4.121, p £ [1, oo], and a > 2. Comparing 



the estimate above with Table [T] we observe that BQCE converges more quickly 
than ATM when a < 3, at the same rate when a = 3, and more slowly when 
a > 3, which we summarize in Table [3j Our numerical experiments concern the 
case a = 3, and we do observe that BQCE and ATM converge at the same rate. 
However, for optimally chosen parameters, BQCE is still much more accurate than 
ATM by a significant constant factor. 



5. Numerical Example 
5.1. Setup of the atomistic model. In the 2D triangular lattice defined by 

*■-(; 

we choose a hexagonal domain J7 a with embedded microcrack as described in Fig- 
ure [T] The sidelength of the domain is TV = 500 atomic spacings, and a de- 
fect is introduced by removing a segment of atoms a t the center of 17 a , A^f := 
{— 5ei, . . . ,5ei} in the microcrack example in Section 5.3 and A!^ := {0, ei} in 



the di-vacancy example in Section 5.4 where the subscripts refer to the number of 
atoms removed in the simulation. 

We supply the domain with Dirichlet boundary conditions by setting y a (0 = 
for £ G A \ f2 a , where F is a macroscopic strain that we will specify below. 
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TOO -50 



50 100 




10 15 20 25 



Figure 1. Deformed configuration in atomic units of a BQCE 
solution for a microcrack in a computational domain with approx- 
imately 3N 2 atoms (N = 100 in this figure, but N = 500 in the 



benchmarks described in Sections 5.3 and 5.4). The color and size 



of the atom positions indicate the value of the blending function. 



The site energy is given by the EAM toy-model (2.2), with 

0(r) = [ e - 2a{r -^ - 2e- a(r -^] , p(r) 
G(p) = C [(p-po) 2 + (p-A)) 4 ], 



- br 



and 



where a, b, c, po £ K are parameters of the model. In all our computational experi- 
ment we choose 



4.4, 



3, 



5, po = 6e~ 



1. 



2.5. 



The model mimics qualitative properties of practical EAM models such as those 



described in 17 34 



To simplify the implementation of the model, we have introduced a cut-off in the 
reference configuration instead of the deformed configuration. In (2.2), we compute 
only sums over first, second and third neighbors. Since there are no rearrangements 
of atom connectivity in our experiments, we do not expect that this affects the 
outcome of the results. 



5.2. Setup of the BQCE method. Our implementation of the BQCE method is 



based on (4.3 1, using standard finite element assembly techniques. The construction 



of the blending function is governed by two approximation parameters: 

• Ko € N denotes the number of atomic layers surrounding the microcrack 
where (3 = 1; 

• Ki £ N denotes the number of atomic layers in the blending region. 



For the microcrack problem we expect a = 3 in the context of Section 4.4 Hence, 
according to Table [I] the choice K\ :— Kq (so that the number of atoms in the 
atomistic region and blending region are comparable) results in the optimal rate of 
convergence for both p = 2 and p = oo . 
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Let <ihop(£, f]) denote the hopping distance in the triangular lattice (with natural 
extension to sets), then we define 

A a := G fl a : d h op(£, A crack ) <K }, and 

A b := {£ G fl a : K Q < d hop (H, A crack ) < K + K 1 }. 

We consider three choices of the blending function, which are all easily defined for 
general interface geometries: 

• QCE: choosing (3 to be the characteristic function of the atomistic region 



A a , and K x = 0, yields the QCE method defined in (4.1). 

• Linear Blending: Let denote the hopping distance from the atomistic 
region, then we choose 

Ai„(£) :=max(l,d(£)/^i). 

• Smooth Blending: The three nearest-neighbour lattice vectors are given by 

ai := ( ) ' ° 2 := ( n/3/2 ) ' and ° 3 := ( rf/l 
Define A?j9(£) := 0(£ + m) - 2/3(0 + /3(£ - a*), and 

Then we define 

ftmooth := argmin{$0S) : = in A a and = 1 in O a \ A a U A b }. 

Roughly speaking, this choice generalizes the cubic spline blending sug- 
gested in |32| . It provides a practical variant for interface geometries that 
are not circular. 

The third approximation parameter is the finite element mesh in the continuum 
region. We coarsen the finite element mesh away from the boundary of the blending 
region according to the rule suggested by the complexity estimates in Table [l] As 
a matter of fact it turns out that the mesh size growth is too rapid to create shape- 
regular meshes, hence we also impose the restriction that neighboring clement can 
at most grow by a prescribed factor; this introduces an additional logarithmic factor 



in the complexity estimates 23 Sec. 7.1 



The resulting energy functional is minimized using a preconditioned Polak- 



Ribiere conjugate gradient algorithm described in 31 . We removed the termination 
criterion in this algorithm and allowed it to converge to its maximal precision, that 
is, until the numerically computed descent direction ceases to be an actual descent 
direction for the energy; this occurs at an accuracy of around 10~ 5 in atomic units. 
We then perform three Newton steps, after which the residual is in the range of 
machine precision, to ensure that the results are not affected by the accuracy of the 
nonlinear solver. 

5.3. Microcrack. In the microcrack experiment, we remove a long segment of 
atoms, Afi f — {— 5ei, . . . , 5ei} from the computational domain. The body is then 
loaded in mixed mode I & I, by setting 
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Figure 2. Convergence rates in the energy- norm (the H 1 - 
scminorm) for the microcrack benchmark problem described in 
Section |5.3| This corresponds to the case p = 2 and a = 3 in 
Table [T] In the above, y & denotes the minimizer of the atomistic 
energy, and j/ ac denotes the corresponding minimizer of the ATM, 
QCE or BQCE problems. 

where F Q oc / minimizes the Cauchy-Born stored energy function W, 7i = 7i = 0.03 
(3% shear and 3% tensile stretch). 

We solve the BQCE problem (to be precise, the reduced atomistic problem 
(ATM) for increasing domain sizes, QCE and BQCE problems for linear and smooth 
blending functions) for increasing parameters Kq, and compute the error relative 
to the exact atomistic solution. 

The relative errors in the W^ 1,2 -seminorm are displayed in Figure [2] the relative 
errors in the W^'^-seminorm are displayed in Figure [3] the errors in the energy are 
displayed in Figure [4] We observe a significant pre-asymptotic regime where the 
rate of convergence is faster than predicted in our theory. For larger computations, 
the rate of convergence approaches closely the predicted rate. 

In particular, it is worth noting that the advantage of a smooth blending function 
only becomes significant in the asymptotic regime and is less pronounced than our 
theory might suggest. Since the precomputation of /3 smoo th is computationally 
cheap and straightforward, we nevertheless propose the choice /9 S mooth for most 
simulations. 

Finally, we remark that the initial dip in the energy error in Figure [4] is most 
likely due to a change of sign in the energy error. 

5.4. Di- vacancy. We perform a second numerical experiment, in order to demon- 
strate that the effects we observed in the microcrack example, which are not covered 
by our theory, are indeed caused by a pre-asymptotic regime. To this end we con- 
sider a di- vacancy defect, where only two neighboring sites A^ 1 := {0, e\} are 
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Figure 3. Convergence rates in the W^' 00 -seminoma for the 
microcrack benchmark problem described in Section [5. 3 1 This cor- 
responds to the case p — oo and a — 3 in Table [I] In the above, 
y a denotes the minimizer of the atomistic energy, and y ac denotes 
the corresponding minimizer of the ATM, QCE or BQCE prob- 
lems. 
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Figure 4. Convergence rates in the energy, for the microcrack 
benchmark problem described in Section |5.3| In the above, E a 
denotes the minimum of the atomistic energy, and E &c denotes the 
minimum of either the ATM, QCE or BQCE energies. 
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Figure 5. Convergence rates in the energy- norm (the H 1 - 
seminorm) for the di- vacancy benchmark problem described in Sec- 
tion |5.4| This corresponds to the case p — 2 and a = 3 in Table [T] 
In the above, y a denotes the minimizer of the atomistic energy, 
and y ac denotes the corresponding minimizer of either the ATM, 
BQCE or QCE problems. 



removed from the lattice. We apply 3% isotropic stretch and 3% shear loading, by 
setting 

where Fq minimizes W , s = 7j = 0.03. 

The relative errors in the W^' 2 -seminorm are displayed in Figure [H] the relative 
errors in the W^'^-seminorm are displayed in Figure |6| the errors in the energy are 
displayed in Figure For the W 1 ' 2 and VF 1,00 -seminorms we now observe precisely 

l— ' 1/9 1 

the rates of convergence DoF ' and DoF predicted by our theory. However, we 
observe a convergence rate DoF -3 / 2 for the energy, rather than the square of the 
rate for the W" 1,2 -seminorm which is DoF -1 . We can, at present offer no rigorous 
explanation for this improved convergence rate, but speculate that it is caused by 
a cancellation effect that is not captured by our theory. 

It is important to point out that, in this second numerical experiment, BQCE 
is not substantially more efficient than ATM. In fact, in the W /1 '°°-seminorm the 
errors are almost identical. This example clearly demonstrates that an atomistic 
coarse-graining scheme need not always lead to improved computational results, but 
that the improvement over a judiciously performed atomistic simulation is highly 
problem dependent. We also remark that for similar examples (e.g., a single vacancy 
with F = F ), we observed that ATM has a higher accuracy than BQCE in both 
the W 1,2 and TT^ 1,00 -seminomas. 
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Figure 6. Convergence rates in the W /1,00 -seminorm for the 
di- vacancy benchmark problem described in Section [5. 4| This cor- 
responds to the case p = oo and a = 3 in Table [I] In the above, 
y a denotes the minimizer of the atomistic energy, and y ac denotes 
the corresponding minimizer of the ATM, BQCE or QCE prob- 
lems. 
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Figure 7. Convergence rates in the energy, for the di- vacancy 
benchmark problem described in Section |5.4| In the above, E a 
denotes the minimum of the atomistic energy, and E ac denotes the 
corresponding minimum of either the ATM, QCE or BQCE ener- 
gies. 
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6. Conclusion 

We have formulated an atomistic-to-continuum coupling, which we call the 
blended energy-based quasicontinuum method (BQCE). Our formulation requires 
only two approximation parameters for the implementation of BQCE: the blending 
function (3 and the finite clement mesh T . For the problems of a microcrack and 



a di-vacancy in a two-dimensional crystal, we utilized theoretical results from 24 
to obtain optimal choices of approximation parameters (blending function and fi- 
nite element grid) to minimize a global error norm and confirmed our analytical 
predictions in numerical tests. 

An interesting open question is how well different atomistic-to-continuum cou- 
plings (in particular, QCE, BQCE and the patch test consistent methods) perform 
if the quantity of interest is localized in the atomistic region. In [27] , this question 
was investigated through a series of numerical experiments in which the Arlequin 
method was used to simulate a one-dimensional chain with a defect. The authors 
concluded that the local error of the Arlequin method could be controlled by moving 
the blending region away from the defect and increasing its size. For the quasicon- 
tinuum method, local quantities of interest were investigated in [2||3j ( 20 . We hope 
to develop rigorous error estimates and optimal approximation parameters for lo- 
calized quantities of interest in a future work. 

At present BQCE is a competitive choice for numerical simulations of atom- 
istic multi-scale problems as it applies to a much broader class than any known 
patch test consistent method. For example, Shapeev's method applies only to two- 
dimensional, simple lattice crystals with pair interactions. Our BQCE method 
can be applied to challenging and physically important problems featuring three- 
dimensional, multi-lattice crystals and arbitrary many-body interaction potentials. 
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